/* 
 * File:   regrid.cc
 * Author: boyan
 *  
 * Created on March 17, 2009, 1:26 PM
 */

#include <math.h>

#include "errlog.h"


#include <iostream>
#include <stdlib.h>
#include <cmath>
#include <cfloat>
#include <string>

// includes from grid3d/
#include "Mesh_grid3d.h"
#include "grid3d.h"
#include "g3d_mesq.h"
#include "Mesh_surf3d.h"

// includes from Mesquite
#include <Mesquite.hpp>
#include <MeshInterface.hpp>
#include <MeshWriter.hpp>
#include <MsqError.hpp>
#include <ShapeImprovementWrapper.hpp>
#include <LaplacianIQ.hpp>

// includes from regrid/
#include "SurfAlign_grid3d.h"
#include "MeshAligner.h"
#include "MultiSurfAlign_grid3d.h"
#include "AlignDomain.h"
// #include "SurfaceAlignDomain.h"
#include "SmartLaplacianIQ.h"

#include "common.h"
using namespace std;

using Mesquite::Vector3D;
using Mesquite::Mesh;
using Mesquite::MsqError;

int mpi_num_procs = 1, mpi_my_id = 0;


static grid3d_t *get_ingrid(int argc, char** argv);
static grid3d_t *make_hold_all_grid(int argc, char** argv, const grid3d_t *ingrid);
static grid3d_t *make_outgrid(grid3d_t *refgrid, LevelsetCode_array &levels);
extern unsigned parse_patches_infile(const char *fname, unsigned *codes,
        char **names, unsigned max_patches, unsigned max_name_length);

int main(int argc, char** argv)
{
#define WRITE_MESH(name_prefix)                                                 \
    {                                                                           \
        FILE *filename = fopen(name_prefix ".bin", "wb");                       \
        save_grid_bin(outgrid, filename);                                       \
        fclose(filename);                                                       \
        Mesquite::MeshWriter::write_vtk(&sparks, name_prefix ".vtk", err);      \
        align_domain.write_vtk(&surf_sparks, name_prefix ".surf.normals.vtk");  \
        align_domain.write_vtk(&sparks, name_prefix ".codes.vtk", true);        \
        cout << name_prefix " quality  " << std::setw(10) <<                    \
                GridQuality(outgrid, outgrid->nodes) << std::endl;              \
    }

    Mesquite::MsqPrintError err(std::cout);

    // load the original grid and construct the domain
    grid3d_t *ingrid = get_ingrid(argc, argv);

    MultiSurfAlign_grid3d align_domain(ingrid);

    {
        unsigned codes[50], num;
        char _names[50][50], *names[50];
        for(unsigned i = 0; i < 50; i++) names[i] = _names[i];

        num = parse_patches_infile(argv[2], codes, names, 50, 50);
        for(unsigned i = 0; i < num; i++)
            align_domain.add_patch(codes[i], names[i]);
    }

    cout << align_domain.patches_to_string();

    align_domain.compile_patches();

    align_domain.write_vtk("align_domain.vtk");

    // generate the hold-all-grid
    grid3d_t *hag = make_hold_all_grid(argc, argv, ingrid);
    MakeInverse(hag);
    {
        FILE *filename = fopen("original_hag.bin", "wb");
        save_grid_bin(hag, filename);
        fclose(filename);
    }

    // wrap a Mesquite::Mesh around hold-all-grid
    Mesh_grid3d hag_mesh(*hag, true); // owns grid

    // this is the surface aligner: decides which nodes from hold-all-grid
    // go inside, outside and on-the-surface
    MeshAligner surf_aligner(&hag_mesh, &align_domain, err);
    surf_aligner.initialize();
    surf_aligner.choose_vertices();
    surf_aligner.dump_levels<Mesh_grid3d::handle2node > ("levelset");

    LevelsetCode_array &surf_levels = surf_aligner.get_levels();
    // update on_bdry codes of hag
    for(size_t q = 0; (int) q < hag->NN; ++q)
    {
        node_t &nd = hag->nodes[q];
        Mesh::VertexHandle h = Mesh_grid3d::node2handle(q);
        if(surf_levels[h] == ls_onsurf)
        {
            Vector3D pos(nd.coords), clo, nor;
            align_domain.snap_point_to_surface(h, pos);
            align_domain.closest_point(h, pos, clo, nor, err);
        }
    }

    // wrap a Mesquite::Mesh around the surface in hold-all-grid
    Mesh_surf3d hag_surf(*hag, surf_levels, ls_onsurf); // does not own
    // for each curve, align surface grid with it.
    std::vector<MultiSurfAlign_grid3d::SurfaceCurve>::iterator C;
    for(C = align_domain.curves.begin(); C < align_domain.curves.end(); ++C)
    {
        // class SurfaceCurve is actually an AlignDomain descendent
        MeshAligner curve_aligner(&hag_surf, &(*C), err);
        curve_aligner.initialize();
        curve_aligner.choose_vertices();

        LevelsetCode_array &curve_levels = curve_aligner.get_levels();
        LevelsetCode_array::iterator R2D2;
        for(R2D2 = curve_levels.begin(); R2D2 != curve_levels.end(); R2D2++)
        {
            if(level_of(R2D2) == ls_onsurf)
            {
                Mesh::VertexHandle h = handle_of(R2D2);
                int n = Mesh_grid3d::handle2node(h);
                Vector3D pos(hag->nodes[n].coords), clo, nor;
                align_domain.snap_point_to_curve(h, pos);
                align_domain.closest_point(h, pos, clo, nor, err);
            }
        }
        char fname[50];
        sprintf(fname, "ls_curve_%d", (int) (C - align_domain.curves.begin()));
        curve_aligner.dump_levels<Mesh_grid3d::handle2node > (fname);
    }

    surf_aligner.move_vertices();

    {
        FILE *filename = fopen("aligned_hag.bin", "wb");
        save_grid_bin(hag, filename);
        fclose(filename);
        Mesquite::MeshWriter::write_vtk(&hag_surf, "aligned_hag_surf.vtk", err);
        align_domain.write_vtk(&hag_surf, "hag_surf.normals.vtk");
        align_domain.write_vtk(&hag_surf, "hag_surf.codes.vtk", true);
    }

    for(size_t q = 0; (int) q < hag->NN; ++q)
    {
        node_t &nd = hag->nodes[q];
        Mesh::VertexHandle h = Mesh_grid3d::node2handle(q);
        align_domain.get_codes(&h, &nd.on_bdry, 1);
    }
    grid3d_t *outgrid = make_outgrid(hag, surf_aligner.get_levels());
    Mesh_grid3d sparks(*outgrid, true); // owns grid

    // from now on we will use align domain with outgrid, not hag anymore
    // let align_domain forget about hag's handles
    align_domain.clear_handles();
    // and tell align_domain about outgrid's handles
    for(size_t q = 0; q < (size_t) outgrid->NN; ++q)
    {
        Mesh::VertexHandle h = Mesh_grid3d::node2handle(q);
        unsigned &code = outgrid->nodes[q].on_bdry;
        if(code)
            align_domain.set_codes(&h, &code, 1);
    }
    for(size_t q = 0; q < (size_t) outgrid->NF; ++q)
    {
        face_t *f = outgrid->faces + q;
        unsigned f_on_bdry = (unsigned) - 1;
        for(int w = 0; w < 3; w++)
            f_on_bdry &= outgrid->nodes[f->nodes[w]].on_bdry;
        Mesh::EntityHandle h = Mesh_grid3d::elem2handle(q);
        if(f_on_bdry)
            align_domain.set_codes(&h, &f_on_bdry, 1);
    }

    // Now I want to improve quality of surface grid.
    // Make a surface Mesh from outgrid
    Mesh_surf3d surf_sparks(*outgrid);
    surf_sparks.set_fixed_onbdry(true);
    surf_sparks.set_fixed_DIR(false);


    WRITE_MESH("outgrid");

    // unmark surface nodes that aren't on curves from being fixed
    std::vector<Mesh::VertexHandle> surf_sparks_handles;
    surf_sparks.get_all_vertices(surf_sparks_handles, err);
    for(size_t q = 0; q < surf_sparks_handles.size(); ++q)
    {
        Mesh::VertexHandle h = surf_sparks_handles[q];
        //if (align_domain.is_handle_on_curve(h))
        //   continue;
        int n = Mesh_grid3d::handle2node(h);
        node_t &node = surf_sparks.grid3d()->nodes[n];
        node.on_bdry = 0;
    }
    // improve quality
    //    Mesquite::InstructionQueue *Queue = new Mesquite::InstructionQueue;
    //    Mesquite::TerminationCriterion *Term = new Mesquite::TerminationCriterion;
    //    Mesquite::LaplacianSmoother *Smoother = new Mesquite::LaplacianSmoother(err);
    //    Term->add_criterion_type_with_int(Mesquite::TerminationCriterion::NUMBER_OF_ITERATES,10,err);
    //    Smoother->set_outer_termination_criterion(Term);
    //    Queue->set_master_quality_improver(Smoother, err);
    //    Queue->run_instructions(&surf_sparks, &align_domain, err);
    //    delete Term;
    //    delete Smoother;
    //    delete Queue;

    SmartLaplacianIQ lap;
    lap.disable_printing_results();
    lap.run_instructions(&surf_sparks, &align_domain, err);

    //    Mesquite::ShapeImprovementWrapper siw(err);
    //    siw.run_instructions(&surf_sparks, &align_domain, err);


    // restore marks of all surface nodes
    for(size_t q = 0; q < surf_sparks_handles.size(); ++q)
    {
        Mesh::VertexHandle h = surf_sparks_handles[q];
        int n = Mesh_grid3d::handle2node(h);
        node_t &node = surf_sparks.grid3d()->nodes[n];
        align_domain.get_codes(&h, &node.on_bdry, 1);
    }
    lap.run_instructions(&sparks, err);

    WRITE_MESH("laped");

    mesq_handle h = mesq_init(outgrid);
    mesq_run_no_wrapper(h);
    mesq_done(h);

    // siw.run_instructions(&sparks, err);

    WRITE_MESH("mesqed");

    // hag and outgrid are released by their owners (hag_mesh and sparks)
    if(ingrid)
        ReleaseGrid(&ingrid);

    return (EXIT_SUCCESS);

#undef WRITE_MESH
}

/****************************************************************************/

/****************************************************************************/

static grid3d_t *get_ingrid(int argc, char** argv)
{
    grid3d_t *grid = NULL;

    if(argc < 2)
    {
        std::cerr << "Missing input grid" << std::endl;
        exit(EXIT_FAILURE);
    }

    string infname = argv[1];

    cout << "Reading ingrid from file " << infname << endl;

    size_t foo = infname.rfind('.');
    if((foo != string::npos) && (infname.compare(foo, string::npos, ".g3b") == 0))
        grid = LoadGridBin(infname.c_str());
        // cout << "extension is " << "g3b" << endl;
    else if((foo != string::npos) && (infname.compare(foo, string::npos, ".g3t") == 0))
        grid = LoadGrid(infname.c_str());
        // cout << "extension is " << "g3t" << endl;
    else
    {
        grid = (grid3d_t *) calloc(1, sizeof (grid3d_t));
        FILE *filename = fopen(infname.c_str(), "rb");
        load_grid_bin(grid, filename);
        fclose(filename);
    }
    if(!grid)
    {

        cerr << "Could not load grid from file" << infname << endl;
        exit(EXIT_FAILURE);
    }

    return grid;
}

static grid3d_t *make_hold_all_grid(int argc, char** argv,
        const grid3d_t *ingrid)
{
    double cmin[3], cmax[3];
    int nints[3];

    double in_h = GridSize(ingrid, ingrid->nodes);
    double h = -1.0;

    string hstr = "-h";
    for(int i = 1; i < argc - 1; i++)
        if(hstr.compare(argv[i]) == 0)
            h = atof(argv[i + 1]);
    if(h < 0.0)
    {
        cerr << "No value of h specified on command line" << endl;
        exit(EXIT_FAILURE);
    }

    node_t *nd = ingrid->nodes;
    for(int j = 0; j < 3; j++) cmin[j] = cmax[j] = nd->coords[j];
    for(int i = 1; i < ingrid->NN; i++)
    {
        nd = ingrid->nodes + i;
        for(int j = 0; j < 3; j++)
        {
            cmin[j] = fmin(cmin[j], nd->coords[j]);
            cmax[j] = fmax(cmax[j], nd->coords[j]);
        }
    }

    for(int j = 0; j < 3; j++)
    {
        cmin[j] -= 3 * h;
        cmax[j] += 3 * h;
        nints[j] = ceil((cmax[j] - cmin[j]) / h);
        if(nints[j] % 2 == 1)
            ++nints[j];
    }

    grid3d_t *grid = GenerateGrid1(cmin[0], cmax[0], nints[0],
            cmin[1], cmax[1], nints[1], cmin[2], cmax[2], nints[2]);

    cout << "MinX = " << cmin[0] << "  MaxX = " << cmax[0] << "  NX = " << nints[0] << endl;
    cout << "MinY = " << cmin[1] << "  MaxY = " << cmax[1] << "  NY = " << nints[1] << endl;
    cout << "MinZ = " << cmin[2] << "  MaxZ = " << cmax[2] << "  NZ = " << nints[2] << endl;

    h = GridSize(grid, grid->nodes);

    cout << "GridSize of ingrid is " << in_h << std::endl;
    cout << "GridSize of hagrid is " << h << std::endl;

    return grid;
}

static grid3d_t *make_outgrid(grid3d_t *refgrid, LevelsetCode_array &levels)
{
    grid3d_t *outgrid = (grid3d_t *) calloc(1, sizeof (grid3d_t));
    int &NN = outgrid->NN;
    // count the number of nodes
    for(size_t nd = NN = 0; nd < levels.size(); ++nd)
        if(levels[Mesh_grid3d::node2handle(nd)] != ls_outside)
            ++NN;
    // allocate nodes
    outgrid->nodes = (node_t*) calloc(NN, sizeof (node_t));
    int *map = (int*) calloc(refgrid->NN, sizeof (int));
    for(size_t nd = NN = 0; nd < levels.size(); ++nd)
    {
        if(levels[Mesh_grid3d::node2handle(nd)] != ls_outside)
        {
            map[nd] = NN;
            outgrid->nodes[NN] = refgrid->nodes[nd];
            ++NN;
        }
        else
            map[nd] = -1;
    }
    // now count the inside elements;
    int &NEL = outgrid->NEL;
    for(int el = NEL = 0; el < refgrid->NEL; ++el)
    {
        elem_t *T = refgrid->elems + el;
        int j, in, out, on;
        for(j = in = out = on = 0; j < 4; j++)
        {
            size_t T_nodes_j = T->nodes[j];
            switch(levels[Mesh_grid3d::node2handle(T_nodes_j)])
            {
                case ls_outside: out++;
                    break;
                case ls_inside: in++;
                    break;
                case ls_onsurf: on++;
                    break;
            }
        }
        // can't have in and out in the same element; can't have all on either
        assert((out * in == 0) && (on < 4));
        if(in > 0)
            ++NEL;
    }
    // allocate elems
    outgrid->elems = (elem_t*) calloc(NEL, sizeof (elem_t));
    for(int el = NEL = 0; el < refgrid->NEL; ++el)
    {
        elem_t &T = refgrid->elems[el];
        for(int j = 0; j < 4; j++)
        {
            size_t T_nodes_j = T.nodes[j];
            if(levels[Mesh_grid3d::node2handle(T_nodes_j)] == ls_inside)
            {
                for(int q = 0; q < 4; q++)
                {
                    assert(map[T.nodes[q]] >= 0);
                    outgrid->elems[NEL].nodes[q] = map[T.nodes[q]];
                }
                ++NEL;
                break;
            }
        }
    }
    free(map);
    EnrichGrid(outgrid);
    UpdateMeasures(outgrid);

    return outgrid;

}






